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Changes in parameters of a physical device can eventually 
give way to catastrophic failure. This paper discusses a pa¬ 
rameter estimation method based on synchronization between 
a model and time series data. In particular, we examine the 
robustness of the method to additive noise in the data. 


I. INTRODUCTION 

Scientists and engineers are often faced with the task 
of estimating the parameters of a physical device using 
time series data obtained from the device, and a math¬ 
ematical model for the dynamics of the device. Because 
the dynamics of a model is governed by its parameters, 
parameter estimation amounts to adjusting the param¬ 
eters in the model until its dynamics accurately mimics 
those of the device. Then, the parameter values of the 
device are assumed to be identical to those of the model. 

Parameter estimation is motivated, in part, by the fact 
that the dynamics of all physical devices change over 
time. If the change arises from normal “wear and tear” 
then it is typically a slow change which eventually gives 
way to a rapid catastrophic failure. Parameter estimation 
aids in avoiding catastrophic failure in two ways. First, 
one can determine regions of parameter space associated 
with failure by estimating and tracking parameter values 
as devices are driven to failure. Second, by estimating 
parameter values of a device at regular intervals (and 
tracking these values) one can determine if and/or when 
the device approaches failure. 

Currently, most devices in the field operate in regimes 
where their dynamics are regular (i.e. periodic or quasi- 
periodic). For these situations, slow changes in param¬ 
eter values can be detected by observing shifts in the 
frequencies and amplitudes of Fourier spectra from time 
series measurements. Recently, there has been specula¬ 
tion about devices whose normal operating regimes have 
complicated and/or chaotic dynamics. Examples of such 
devices are bearings under high stress, electronic compo¬ 
nents with nonlinear response functions, mixing and elec¬ 
trochemical reactions, and high speed cutting tools. For 
these examples, complicated dynamics are unavoidable, 
and may even be desirable. More importantly, Fourier 
spectra will not exhibit sharp peaks. Therefore, other 


methods for estimating and tracking changes in, param¬ 
eter values are needed. 

This paper discusses parameter estimation for sys¬ 
tems whose underlying dynamics are complicated and/or 
chaotic. An early example of parameter estimation in 
the nonlinear dynamics literature is Hudson et.al. 0. In 
this work the authors constructed neural net models that 
mimicked a few of the period doubling bifurcations found 
in an electrochemical experiment. 

Recently, Baker et.al. proposed a least squares ap¬ 
proach that requires time series vectors which are simul¬ 
taneous measurements of all of the phase space variables. 
The approach also requires numerical approximation of 
derivatives from the time series. Horbelt et.al. (build¬ 
ing on work by Baake et.al. Q) implemented a shooting 
method. They used the method, and experimentally ob¬ 
tained data, to estimate parameters for a model of cal¬ 
cium release in muscle cells . Also, a novel method using 
symbolic dynamics has been developed for systems that 
exhibit complicated temporal or spatio-temporal dynam¬ 
ics (§. 

The approach discussed in this paper estimates the pa¬ 
rameters of a physical device by selecting values which 
yield “the best synchronization” between a time series 
from the device and a mathematical model for the de¬ 
vice. The method is related to Kalman filtering and 
relies on two ideas. First, if two dynamical systems are 
identical then it is possible for them to synchronize (i.e., 
both systems follow identical phase space trajectories). 
Second, if this behavior is stable then (for the noise free 
case) deviations from synchronous behavior converge to 
zero as differences in parameter values between the two 
systems converge to zero. These ideas imply that (in the 
absence of modeling errors) the best estimate for the pa¬ 
rameter values of a device are those of the model which 
best synchronizes to a time series from the device. 

The synchronization approach to mrameter estimation 
has been discussed by Parlitz et. al. [Qlg as well as Mayb- 
hate and Amritkar . The work discussed in this paper 
focuses on how noise in experimental data effects the ac¬ 
curacy of parameter estimation. Refs. [M either did 
not discuss noise, or only examined smal noise levels. 

Furthermore, our implementation differs from the ones 
in previous papers. For example. Refs. [^,|[P used vari- 



ous down hill search methods to minimize the cost func¬ 
tion. However, it is known |5|,[l^ (and we have observed 
in the numerical experiments discussed below) that a cost 
function can have many local minima where simple down 
hill methods can be trapped. We avoid this problem in 
some of our work by using annealing to search the pa¬ 
rameter space. 

Another important difference between our approach 
and Refs. @] and concerns the time series used in 
the estimation. The previous authors assume the exper¬ 
imental measurements are continuously recorded. Thus, 
the time series has measurements at every value of time. 
When performing numerical experiments with this type 
of time series it is reasonable to replace the time series 
by a set of ODE’s. 

In our work we assume the time series consists of exper¬ 
imental measurements obtained at discrete increments of 
time (given by a sampling interval). When performing 
numerical experiments with this type of time series one 
can not replace the time series by a set of ODE’s. Rather, 
one must use an approach that explicitly accounts for the 
discreteness of the time series. 

Furthermore, Ref. assumes that one of the phase 
space variables can be measured. We only assume a 
scalar time series corresponding to some experimental 
measurement. It does not have to be one of the phase 
space variables. Finally, our numerical experiments si¬ 
multaneously estimate several parameters, while Ref. ^ 
typically estimates only one or two parameters for a three 
dimensional model. We have tested our methods on a 
high dimensional systems ( an 11 dimensional general¬ 
ized Rossler model). 

Although modeling errors are inevitable in any math¬ 
ematical description of a device, this paper does not 
address this issues. We conjecture that noise is the 
dominant obstacle associated with parameter estimation. 
Our examples represent electronic circuits, chemical re¬ 
actions, and a high dimensional version of the Rossler 
system. The results indicate that this approach to pa¬ 
rameter estimation holds promise for real applications. 

The remainder of this paper is organized as follows. In 
Section we define and discuss synchronization between 
two or more dynamical systems. In Section III we discuss 
parameter estimation via synchronization. Examples and 
the effects of noise are given in Section [V. Conclusions 
are presented in Section 0. 


II. SYNCHRONIZATION 

In this section we define and discuss synchronization 
between two coupled deterministic dynamical systems. 
The numerical experiments use two different synchroniza¬ 
tion techniques, each of which is discussed. The func¬ 
tional forms of the uncoupled systems are identical, so 
they only differ in the values of their respective param¬ 
eters. The coupling between the systems is denoted by 


E{x,y), where £c,y S 1R‘^ are the instantaneous states 
of the two systems. The equations of motion for this 
arrangement are 


dx 

dt 

dt 


F{x;p*) 

F{y;p) + E{x,y). 


( 1 ) 

( 2 ) 


Equation ( 0 ), the drive system, represents the physical 
device, while Eq. (|^, the response system, represents a 
mathematical model of the device. The parameters of 
the device and model are represented by p*,p G IR", 
respectively. The coupling shown in Eqs. ( 0 ) and (|^) 
is called drive-response coupling, and is natural for our 
application. If the coupling in Eq. (^ vanishes when 
X — y (i.e., E{x,x) = 0) then synchronization is defined 
as a;(t*) = y(t*) at some time —oo < < oo. If this 

occurs then, in the absence of noise, determinism implies 
that x{t) = y{t) for t> t^. 

However, the systems we examine are chaotic, no 
model of a physical device can ever be exact, the state of a 
model can never exactly match the state of a device, and 
noise is never absent from real data. Therefore, synchro¬ 
nization between uncoupled chaotic systems never really 
exists because approximation synchronization {x ~ y) 
between chaotic systems is always unstable. 

These remarks imply that the coupling between 
Eqs. ( 0 ) and ( 0 ) is necessary to stabilize synchronous 
motion. Over the last several years it has become in¬ 
creasingly clear that ii p = p* then there are many 
choices of E{x,y) which will result in stable synchro¬ 
nization 1^. For these types of coupling, if ||a; — y[| is 
small then limt^oo \\x — y\\ — 0, where || • || denotes the 
Euclidean norm. 

One form of coupling that satisfies E{x,x) = 0 and 
often yields stable synchronization is 


E{x, y) = K[h{x) - h{y)], (3) 

where K (the so called gain vector) could have only one 
nonzero component. Thus, the coupling between the 
drive and response systems could involve only one com¬ 
ponent of F{y;p). 

It is usually impossible to simultaneously measure all 
of the state variables of a physical device. Typically, one 
experimentally measures a scalar function of the state of 
the device. To make contact with this more realistic case 
we introduce the measurement functions, h{x) and h{y). 
These functions represent the experimental measurement 
process which transforms the instantaneous state of the 
device or model into a scalar measurement. 

The time series, h[x{tn)], represents scalar data mea¬ 
sured at discrete times = nAt, n = 1, 2,..., where At 
is the sampling interval. In contrast, numerical integra¬ 
tion of the model requires E at arbitrary values of t. To 
overcome this problem we have used two different meth¬ 
ods for calculating E at times not equal to one of the 
tnS. One method selects all data within a time window 




centered at t. Least squares is used to fit the data points 
within this window to a polynomial, and the polynomial 
is used to approximate h[x{t)] [Q. For this method the 
response system, Eq. (||), is continuously forced by the 
time series. 

The second method is called sporadic driving [^ . 
Here, E{x,y) = 0 if t is not equal to one of the t„’s. 
If t is equal to one of the t„’s then E is given by Eq. i)- 
Eor this method, the response system (the model) is dis- 
continuously kicked by the time series at multiples of the 
sampling interval, and evolves freely between kicks. 

Although the two methods differ in their implementa¬ 
tion, each leads to stable synchronization. Furthermore, 
each is a natural choice when confronted with a drive 
system that is a time series instead of ODE’s. 

III. PARAMETER ESTIMATION 

Assume we are given a time series, h[x{tn)] for n = 
1, 2,N, representing experimentally measured data. 
Our task is to select values for the model parameters such 
that its dynamics are the same as those of the device 
which generated the data. If the parameter values of the 
device are p*, and the estimated values are p, then the 
goal of parameter estimation is to achieve p — p* ■ 

The synchronization method uses coupling that guar¬ 
antees stable synchronization when p = p*, and the drive 
system is an ODE. This amounts to determining K in 
Eq. (^. To hnd p we choose parameters, p, which mini¬ 
mize 


time series, and calculated the distance between it and 
its nearest neighbor in the attractor reconstructed from 
the experimental data. If the dynamics of the model and 
the device are the same then every point on the attrac¬ 
tor reconstructed from the model time series will have 
a near neighbor on the attractor reconstructed from the 
experimental data. More importantly, the distance to 
the nearest neighbor will not be large. 

However, if the dynamics of the model and the device 
are different then some points on the model attractor will 
be far away from all points on the experimental attractor, 
and the distance to the nearest neighbor will be large. If 
this is the case then we choose a new initial condition in 
parameter space and minimized Eq. (i again. 

The last example used a gradient descent method 
(POWELL in Ref. |^) to search the parameter space [||. 
To avoid a local minima, we change each parameter by 
an amount of order ~ 10“^, and the descent is started 
again. This procedure is repeated until we find a value 
for p that is stable against such perturbations. 

Numerical experiments indicate that we typically find 
a value of p that is near p* using only one or two initial 
conditions in parameter space. (The annealing proce¬ 
dure we use is probably not the best implementation of 
annealing. We believe that a more sophisticated imple¬ 
mentation would remove this problem p(^ .) For the pas¬ 
sivation example discussed below, the gradient descent 
method frequently fails to find the global minimum on 
the first try. 

The numerical experiments use the following coupling 
between the drive and response systems. 


x{p,p*) 


1 

— ^(/i[a;(n)] - h[y{n)]f 


( 4 ) 


Equation (^ explicitly notes that deviations between 
measured output from the model and measured output 
from the devices are functions of the parameters of the 
model, p, and the device, p*. This dependence arises 
from the dependence of the trajectories, x and y, on 
the parameters in Eqs. (|^) and (||). We denote the 
value of p which minimizes Eq. (jj) by p, and note that 
x{p*tP*) = 0 in the absence of noise. 

The numerical experiments typically use M/N ~ 1/3 
or 2/3. The remaining N — M points are treated as a 
transient phase used to initialize synchronization. 

The first two examples used annealing to search param¬ 
eter space and avoid becoming trapped in local minimum 
far from p* ||l^ . If p is from such a local minimum then 
we expect the dynamics of the model to be very different 
from that of the device. To detect this, we performed a 
simple post-minimization test. After finding p we gener¬ 
ated a time series, h{y), from the model and used time 
delay embedding of h{x) and h(y) to reconstruct attrac¬ 
tors for the data and the model ||^,|^ . 

To test the similarity of the attractors, we selected 
a point on the attractor reconstructed from the model 


E{x,y) = K[{h{x)+r]) - h{y)]. (5) 

Here, y represents noise, and {h{x) + rf) represents an 
experimental measurement that has been contaminated 
with additive noise. In the numerical experiments y = 
ctI]A( 0, 1), where N{Q, 1) denotes random numbers nor¬ 
mally distributed with mean zero and unit standard de¬ 
viation. E represents the size of the experimental signal, 
and was equated to the standard deviation of the clean 
time series, h{x). 

We could have examined dynamic noise. However, dy¬ 
namic noise is equivalent to modeling errors which change 
instantaneously in time. Since we have chosen not to ad¬ 
dress modeling errors it is reasonable to ignore this type 
of noise. 


IV. EXAMPLES 

We now discuss the effect of noise on the accuracy of 
parameter estimation. The first example used a model of 
the electronic circuit shown in Fig. |^. This circuit consist 
of a nonlinear amplifier, V, which transforms input volt¬ 
age, X, into output, af{x) [Q. The parameter, a, char¬ 
acterizes the gain around x = 0. The amplifier has linear 
feedback consisting of a series connection to a low-pass 





FIG. 1. A block diagram for an electronic circuit modeled 
by Eqs. (|) and (|). 


filter, i?C", and a resonant circuit LC. The system has 
been studied experimentally, is known to exhibit chaotic 
dynamics, and has been shown to synchronize with a copy 
of itself for many coupling schemes . 

It is known that a good model for this circuit is the 
following system of ODE’s 

dxi 

^ = -Xi - 6X2 + X3 (6) 

at 

^ = 'y[F{xi) -Xa] - PX 2 
with F{x) defined by 

( 0.528a X < -1.2, 

F(x) = < ax{l-x^) -1.2 < a; <1.2, (7) 

[ -0.528a a; > 1.2 

The numerical experiments use Eqs. d) and (|^) for the 
response system, continuous coupling given by Eq. i), 
and parameter values p* = [^, 7, /3, a] with 6 = 0.43, 
7 = 0.1, [3 = 0.72, and a = 16. For this circuit it is 
straightforward to measure the voltages associated with 
X. However, to make contact with other experimental 
devices we used the following (arbitrarily selected) mea¬ 
surement function 

/i(a;) = 2 xi-I- 3 a: 2 -b 0.5x3. ( 8 ) 

The gain vector is K = [0, 0,10]. In the absence of noise 
this coupling leads to stable synchronization when the 
drive and response systems are ODE’s. The sampling 
interval is At = 0.5 and the data is shown in Fig. |(a). 
For this example, a, the size noise ranged from a low of 
0.001 (representative of very clean data) to a high of 1 
(representative of very noisy data). 

We use an Adams-Bashforth-Moulton code numeri¬ 
cally integrate the model j^. The initial guess for p 
in the annealing procedure is selected at random from 
within the range pi S [0,1], p 2 G [0,2], ps G [0,5] and 
P4 G [0,20]. In the absence of noise, the method cor¬ 
rectly estimated p = p* to within machine precision. 


The numerical experiments with noise examined 100 dif¬ 
ferent value of O'. For each value we calculate p. The 
results of the numerical experiments (see Fig. ^(b)) indi¬ 
cate that lli? — p*|| oc cr||p*|| over much of the range of 
noise strengths. This result confirms earlier theoretical 
predictions [^ ). 

The next example uses a model for metal passivation. 
Passivation is a loss of chemical reactivity associated with 
metal corrosion in aqueous solutions. It is an electro¬ 
chemical reaction that exhibits the full zoo of nonlin¬ 
ear dynamical behavior, including periodic cycles, period 
doubling, multistability, hysteresis, chaos. These behav¬ 
iors are not unusual in electrochemical reactions [pT|-p^ . 

We use the following passivation model 


dt 

dOpH 

dt 

ddp 

dt 


a0M — bY 


Y9M — [r + exp(—/30 oh)]0oh + ^sOpOyi 


xOpn - S^o^M, 


(9) 


where 9m = (1 — 6*0 — ^oh), and the parameters p* = 
[a, b, r, s] are a = 2 x 10“"^, b = 10“^, r = 2 x 10“^, 
and s = 9.7 x 10“®. We fixed the last parameter at 
/? = 5. This example is subtle because of the small size 
of the parameters and the attractor. To overcome this we 
rescale the variables by subtracting the mean and divid¬ 
ing by a constant roughly equal the width of the attractor 
in that coordinate. Time is also rescaled by a factor of 
10,000 [^. The coupling between drive and response 
system uses 


h(x) = 0.2xi-b X2-b 0 . 5 x 3 . (10) 

with gain vector K = [0,100,0]. Once again h{x) is an 
arbitrarily chosen function, and in the absence of noise 
this coupling leads to stable synchronization when the 
drive and response systems are ODE’s. The sampling 
interval of the numerical experiments is At = 0.01 (in 
rescaled time) and the data is shown in Fig. ^(a). 

For this example, cr, the size noise ranged from a low 
of 0.001 (representative of very clean data) to a high of 1 
(representative of very noisy data). 100 different value of 
a were examine. The results, shown in Fig. ^ (b), indi¬ 
cate that the relative accuracy of parameter estimation 
is two digits or more up to cr ~ 0.1 and about 30 % for 
0 dB noise. This surprisingly robust result implies that 
synchronization is worthy of further investigation as a 
parameter estimation technique. 

The final example is a high dimensional generalization 
of the Rossler system [p^ 

dxi 

—rr = -xi -b axi 
dt 

^ = Xj_i - Xj+i (j = 2, ..., L - 1) 

^ = bFdxL(xL-\ -c). 























FIG. 2. The electronic circuit model, (a) Short samples 
of the time series used in the numerical experiments. The 
noisy time series represents 0 dB noise, (b) The relationship 
between the relative error in the parameter estimates and the 
size of the noise. Here, Ap = p — p*. 


FIG. 3. The metal passivation model, (a) Short samples 
of the time series used in the numerical experiments. The 
noisy time series represents 0 dB noise, (b) The relationship 
between the relative error in the parameter estimates and the 
size of the noise. Here, Ap — p — p*. 






















Here, p* = [a, b, c, d\ with a = 0.29, b = 0.1, c = 2, 
and d = 4. This system is the well known Rossler sys¬ 
tem with an additional L — 3 degrees of freedom. It has 
a stable chaotic attractor when L is odd. Our numer¬ 
ical experiments used L = 11, yielding a hyperchaotic 
attractor with a Lyapunovdimension of Hz, ~ 10.1. For 
this example we used for synchronization the method of 
sporadic driving with h{x) = xi as the drive signal and 
a sampling interval of At = 1. The results were calcu¬ 
lated for 15 different noise levels ranging from a = 10“® 
to cr = 0.63. The results are shown in Fig. As in the 
other examples, we observe a linear scaling of the param¬ 
eter relative error with the noise strength a and find 
good estimates up to u ~ 0.15 The quality of the param¬ 
eter estimates in relation to the noise strength is almost 
the same as for the electronic circuit example. This result 
confirms that the method can be applied to highdimen¬ 
sional systems, even if only very noisy time series data is 
available. 


V. CONCLUSION 

In this paper we have examined the effects of noise 
on the accuracy of parameter estimation for physical de¬ 
vices. The parameter estimation method we examined 
uses experimentally measured time series and a mathe¬ 
matical model of the device. The estimates returned by 
the method are those which yields the smallest deviation 
from synchronization between the dynamics of the model 
and the measured time series. 

Numerical experiments indicate that this method is 
robust to additive noise in the time series. We find that 
even at 0 dB signal to noise ratio we can still obtain 
reasonable accuracy in the parameter estimates. The 
method works also for noisy time series coming from high¬ 
dimensional systems. The obvious next step is to test 
this method “at sea”, by using experimental data com¬ 
ing from electronic [||, physical, chemical, ... systems. 
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